Ballistic deposition patterns beneath a growing KPZ interface 
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We consider a (1 + 1) dimensional ballistic deposition process with next-nearest neighbor interac- 
tion, which belongs to the KPZ universality class, and introduce for this discrete model a variational 
formulation similar to that for the randomly forced continuous Burgers equation. This allows to 
identify the characteristic structures in the bulk of a growing aggregate ( "clusters" and "crevices" ) 
with minimizers and shocks in the Burgers turbulence, and to introduce a new kind of equipped 
Airy process for ballistic growth. We dub it the "hairy Airy process" and investigate its statistics 
numerically. We also identify scaling laws that characterize the ballistic deposition patterns in the 
bulk: the law of "thinning" of the forest of clusters with increasing height, the law of transversal 
fluctuations of cluster boundaries, and the size distribution of clusters. The corresponding critical 
exponents are determined exactly based on the analogy with the Burgers turbulence and simple 
scaling considerations. 



PACS numbers: 02.50.-r, 05.10.-a, 05.40.-a 

I. INTRODUCTION 

Over the past few decades, the problem of growth of 
aggregates by sequential stochastic deposition developed 
into one of the most extensively studied topics in statisti- 
cal physics 1]. Much effort has been put into theoretical, 
numerical, and experimental investigation of the result- 
ing patterns. Several theoretical models have been pro- 
posed, including the famous Kardar-Parisi-Zhang (KPZ) 
[2] and Edwards-Wilkinson (EW) Q models, the Re- 
stricted Solid-on-Solid (RSOS) (| and Eden @ mod- 
els, the models of Molecular Beam Epitaxy (MBE) Q, 
Polynuclear Growth (PNG) 04111. 

and several ramifi- 
cations of the Ballistic Deposition (BD) model [12H15] ]. 
Within the latter, in the simplest setting, one assumes 
that elementary units ("particles") follow ballistic tra- 
jectories in space and adhere sequentially to a growing 
aggregate ("heap"). Despite its extremely transparent 
geometric formulation, the problem of stochastic growth 
still is one of the most puzzling problems in statistical 
mechanics. 

The available theoretical analysis of stochastic deposi- 
tion focuses almost exclusively on the enveloping surface 
h(x, t), involving a statistical study of its height distribu- 
tion and the corresponding scaling exponents. Here we 
quote just a few prominent results. The essential scaling 
relations characterizing the growing aggregate are 

(Varh{x,t)) 1/2 ~ t 1/3 , 
(h(x,t)h(x + t 2 / 3 l,t))-(h) 2 ~ t 2 / 3 F{l). 



Here h(x, t) — h(x, t) — ct, c — limt_>oo t^ 1 h(x, t) is an av- 
erage speed of growth, and F{1) is a rescaled correlation 
function. The exponents 1/3 and 2/3 were determined 
already in Q for the KPZ model and then observed in 
a variety of other growth models. Then in Refs [1, [l6| 
it was realized that the distribution of a rescaled PNG 
height i _1 / 3 (/i(0, t) — 2t) converges as t — > oo to the 
Tracy- Widom distribution [18| for the Gaussian unitary 
ensemble (GUE), which appears in the theory of ran- 
dom matrices. Moreover, the full rescaled PNG surface 
t" 1 / 3 (/i(a;t 2 / 3 , t) — 2t) + x 2 converges to a version of the 
Airy stochastic process Airy 2 (x) [9\ whose one-point dis- 
tributions are precisely Tracy-Widom. Distribution of 
maximal heights of the (1 + 1) dimensional Edwards- 
Wilkinson and KPZ interfaces has been determined ex- 
actly in Ref. [r7| . 

It should be pointed out that the BD model considered 
below involves the point-to-line last-passage percolation 
while PNG corresponds to the point-to-point setting [111 ]. 
Correspondingly the limit processes are different: it is 
Airyi for the point-to-line and Airy 2 for the point-to- 
point. Note also that the one-point distribution for the 
Airyi process is given by the Gaussian orthogonal en- 
semble (GOE) distribution rather than the GUE distri- 
bution, which corresponds to Airy 2 . 

Since a similar convergence to Airy processes is ob- 
served in other growth models such as TASEP (28j , it is 
becoming customary to speak of "the KPZ universality 
class" whenever such limit behavior is present. For ex- 
ample, a KPZ scaling has been shown in Refs fToj - f^lj l for 
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FIG. 1: Snapshot of a heap obtained by ballistic deposition 
of N = 2000 particles in a periodic box of size L = 100 
with next-nearest-neighbor interactions. Black lines trace the 
channels ("crevices") between adjacent clusters. 



the BD model in the thermodynamic limit. 

Much less is known, however, about the structure of 
BD patterns beneath the enveloping surface. Here is just 
one puzzle: analytic arguments 22] predict that the ex- 
pected density p sur f of local surface maxima in a (1 + 1) 
dimensional ballistically growing heap is psurf = 1/3, 
whereas extensive numerical simulations show that the 
mean bulk density, pbuik, of the (1 + 1) dimensional heap 
is about 0.25 « 1/4. To date there is no satisfactory 
quantitative explanation of this mismatch. 

The statistics of the growing heap are determined by 
its striking internal structure, revealed in numerical sim- 
ulations as well as in the recent experimental analysis of 
electrochemically formed silver branched patterns [23| . 
This structure consists of a "forest" formed by tree-like 
clusters of different size, which are separated by a dual 
network of tree-like channels or "crevices" (Fig. [T}. 

As the heap grows, clusters randomly collect parti- 
cles and thus spread and isolate their neighbors from the 
"rainfall" of incident particles, suffocating their growth. 
Consequently the number of clusters present at height 
h in a growing aggregate is a decreasing function of h. 
We remark that this "suffocation" mechanism, as well 
as the growth patterns in the BD model, bear certain 
similarity to those observed in diffusion limited aggre- 
gation in a hard-core lattice gas on a zero-temperature 
boundary |24| , although the two models belong to differ- 
ent universality classes and their quantitative behaviors 
are in no direct correspondence. 

In the present work we undertake investigation of clus- 
ters and crevices based on a novel systematic analogy 
with turbulent structures in randomly forced Hamilton- 
Jacobi equations. This allows us to conclude that BD 
belongs to a large group of models within the KPZ uni- 
versality class, such as PNG model, TASEP, and others. 



It turns out that BD like other models mentioned above 
admits a variational formulation. Moreover the analogy 
with Hamilton-Jacobi dynamics enables us to suggest a 
novel concept of equipped Airy process, a buildup on top 
of classical Airy processes which also takes into account 
the geometrical structure of the optimal paths (maximiz- 
ers of the action, see Section HIT]) . The random field of 
optimal paths arises naturally in the context of stochas- 
tically forced Burgers equation 1251 126| . 

In a recent experimental work [23J the size distribution 
P{m) of frozen structures formed by electrochemically 
grown silver branching patterns has been analyzed. The 
authors found that the probability P(m) to have a cluster 
of size m exhibits scale invariance, i.e. P(m) ~ m T , with 
a critical exponent r = 1.37 ± 0.04. In our work we com- 
pute this exponent analytically (r = 7/5) and show that 
the scaling behavior conjectured in |23[ actually holds as 
well as two other power laws governing the "thinning" 
of the forest of clusters with increasing height and the 
transversal fluctuations of the cluster boundaries. 

The paper is organized as follows. In Section [TT] we 
specify the model and define its main structural features. 
Section Mil contains an analysis of the structural simi- 
larity of BD patterns to "minimizers" and "shocks" in 
the Burgers turbulence (26|, based on the common vari- 
ational formulation of the two models. In Section IIVI 
we discuss the KPZ scaling in the BD model and intro- 
duce the notion of an equipped Airy process. Building on 
these developments, in Section [V] we compute the main 
scaling exponents of the BD model. Section IVII contains 
concluding remarks and outlook for future work. 

II. THE MODEL AND BASIC DEFINITIONS 
A. The NNN ballistic deposition model 

A standard (1 + 1) dimensional BD model with next- 
nearest-neighbor (NNN) interactions can be formulated 
as follows (see also Refs [l9l - l2l1 |). Consider a box divided 
into L columns of unit width each, enumerated with in- 
dex i (i = 1,2, ...,L). For simplicity we assume the 
periodic boundary conditions, so that the leftmost and 
the rightmost columns are neighbors, and identify the 
index value with L. 

At the initial time t = the system is empty. Then, 
at each time step t — 1,2,..., i max , an elementary unit 
("particle") of height £ and width 1 is deposited at a 
column i(t) chosen randomly with uniform distribution. 
Define 

Vi(t) = { h * = (la) 

As shown in Fig. [TJ particles deposited in adjacent 
columns interact in such a way that they can only touch 
each other at corners or at top and bottom, but never 
along their vertical sides. Let the height of column i at 
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time t— 1 be hi(t—l). Upon adding a particle it changes 
according to 

(max{/ii_i(t - 1), h t (t - 1), h i+1 (t - 1)} + £, 
%(*) = !, (lb) 
/li(t-l), r?i(t) = 0. 

This dynamics is supplemented with the initial condition 
/i 4 (0) = for all 1 < i < L. Eqs (JTaJ), ^TbJ completely 
describe updating rules for the NNN discrete ballistic de- 
position. 

We will use Eq. (| represented in a different form. 
Define the "thin" and "thick" discrete "5-functions" 

= Joe \k-i\>0, t = joo \k-i\>l, 
w (0 |fc-i| =0, M [0 \k-i\ < 1 ' 

Consider first the trivial dynamics described by the equa- 
tion hi(t) = h%(t — 1). It can be rewritten as ft,j(t) = 
max [/ife(t — 1) — L\ J: indeed, max [/ifc(t — 1) — oo] = — oo 

k k=fci 

and therefore hi(t) = max{hi(t — 1),— oo} = hi(t — 1). 
It is now clear that the stochastic equation (|lb[) can be 
recast in the form 

hi(t) = max [h k (t - 1) - L v k f 3 ] + &h(t). (3) 

This dynamics should be compared with the commonly 
used discrete equation with "additive noise" describing 
the (1 + 1) dimensional polynuclear growth Q, which in 
our notation takes the form 

hi{t) = max [h k (t - 1) - L\ J + Afc(t).l' (4) 

fc 

According to Eq. the height ft,, remains unchanged 
(quenched) if nothing is deposited to column i at time 
t. On the contrary, in Eq.Q the height hi relaxes spon- 
taneously even in the absence of deposition to column i 
at time t because hi(t) is defined to be the maximum of 
the triple {hi-i(t - 1), hi{t - 1), h i+ i(t - 1)}. Note that 
process described by Eq.Q is sometimes referred to as 
"dynamics with multiplicative noise." 



B. Clusters, crevices, and scaling exponents 
in the growing heap 

Let us now take a closer look at Fig. [TJ We say that 
two particles in a heap are connected if they touch one 
another at corners or if one is situated directly on the 
top of the other. 

It often happens that the upper particle is connected 
simultaneously to two lower particles. For reasons that 
will become clear shortly, it is better to avoid these "one- 
on-two" configurations. The model is therefore slightly 
augmented: one assumes in Eqs (jlb[) and ([3]) that 



where are independent normal random variables. It 
is clear, and well supported by numerical experiments, 
that this modification removes the possibility of "one-on- 
two" configurations while preserving, within the limits of 
statistical errors, statistical characteristics of the heap for 
1=1. Alternatively one might resolve "one-on-two" con- 
figurations for £ = by simply disconnecting the upper 
particle from one of its two lower neighbors at random. 
Either way, elimination of one-on-two configurations al- 
low us to define a unique "path" corresponding to every 
particle, namely a backward directed chain of connected 
particles going from a given particle to the bottom level 
of the heap. 

Consider all connected paths originating from the top- 
most particles. These paths can merge. We define the 
backbone of a cluster as the connected set of such paths, 
i.e., the union of all paths that end up at the same bot- 
tom level particle. It is easy to see that the bottom level 
particles are split into two classes: those that are reached 
by the paths originated at the top of the heap and those 
that are not. Obviously the first class gets smaller as t in- 
creases. For every particle from this class define cluster 
as the collection of all paths ending up at this particle. 
The difference between a backbone and a cluster is that 
clusters contain paths not necessarily originating from 
the top level particles. 

We say that a pair of two top level particles occupying 
adjacent columns defines a shock, which is located be- 
tween them, if they belong to two different clusters. The 
channel of white space between two neighboring clusters 
is called a crevice. Clearly every crevice is associated 
with a shock at the top, and the connected paths from 
top particles defining the shock form the left and right 
boundary of a crevice. Connecting shocks at adjacent 
time moments, we get curves that branch forward in time 
and play a role dual to that of backbones. These curves 
are sketched in Fig. [T]in black. 

It is clear from Fig. [T] that many channels that are ini- 
tially present at bottom of the bulk then merge at some 
height, blocking the growth of the clusters situated in be- 
tween. Thus crevices have tree-like structure just as clus- 
ters, but contrary to clusters they merge upward. This 
causes the number of percolating clusters and crevices to 
decrease as a function of h. In the thermodynamic limit 
this behavior is characterized by the following three scal- 
ing exponents whose values are identified in Section fVl 

The thinning exponent a characterizes the expected 
number (c(h)) of percolating crevices (or, equivalently, 
percolating clusters) at height h: 

(c(h)) ~ h~ a (6) 

The wander exponent (3 characterizes the expected 
mean square displacement (in the units of L) of the 
boundary of percolating cluster between the bottom of 
the bulk and a specified height h: 



i = e(t) = i + io- w £{t), (5) 



(Ax 2 (h)) - h p 



(7) 
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The mass exponent r characterizes the mass distri- 
bution of clusters: 

P(m) ~ m- T , (8) 

where P(m) is the proportion of clusters of mass m in 
the ensemble. 



III. BALLISTIC DEPOSITION AND BURGERS 
TURBULENCE 

A. Variational formulation of the BD 

The discrete equation 

hi(t) = max [h k (t - 1) - L kii (t)} + trn{t), (9) 

k 

whose particular cases for specific choices of Lks{t) are 
the BD model © and the discrete PNG model g]), ad- 
mits a natural variational formulation. 

Fix some initial condition hi(0) and consider the 
discrete "variational" problem of finding a trajectory 
(7(0), 7(1), . . . ,7(i)) that satisfies the "boundary condi- 
tion" j(t) = i and maximizes the discrete "action" 

■4o(7) = /i 7 ( )(0)- I L 7(s-i),7(s)( s ) -^7(s)( s )]- 

l<s<t 

( 10 ) 

The function in the square brackets plays a role of a dis- 
crete "Lagrangian" of the system. The problem bears an 
obvious resemblance to the zero temperature limit of the 
free energy of a statistical system, expressed as the sum 
over configurations 7: 

lim T\n(e^ Al H he^ N ) ->■ max-TAi, . . . , A N ] 

Another obvious connection is with mechanics, where the 
dynamical trajectory can be found by optimizing the cor- 
responding action (in our case, at variance with the usual 
convention, the action is maximized). 

Action maximization in Eq. ([10[) is related to solving 
Eq. © as follows. To be specific, consider the BD growth 
© , where particles are added to the system as "dropping 
events" (i(s),s) in (1 + 1) dimensional discrete space- 
time. Maximization of the action A l Q in (|10l) amounts to 
finding a trajectory that terminates at {i,t) and passes 
through a maximal number of dropping events under 
the following constraint: the trajectory stays constant, 
7(s) = 7(s — 1), unless j(s — 1) = i(s) ± 1, i.e., there is a 
dropping event in adjacent column. In the latter case the 
trajectory may (but does not necessarily have to) jump 
to i(s) at time step s. Note that for the PNG model 
(|4]) this constraint is relaxed: a trajectory may jump at 
all times, but only to adjacent columns. Otherwise the 
two models are structurally similar, and the rest of the 
argument in this subsection applies to both. 

The lack of a strict obligation to pass through an adja- 
cent dropping event allows to "collect" dropping events 



more efficiently: it is easy to construct trajectories for 
which it is more profitable, from the point of view of 
maximizing the number of dropping events, to skip some 
isolated dropping events in order not to be driven away 
from a later series of several adjacent dropping events. 

Direct maximization of the action (jlOl) is a difficult 
problem because the solution depends on the whole fu- 
ture history of dropping events. Observe however that 
for all 1 < j < L, 1 < s < t the height function hj(s) 
gives the maximal number of dropping events available 
for a trajectory coming to the point (j,s), and this fact 
can be exploited to construct a maximizing trajectory in 
reverse time. 

Consider again the BD case where Lkj(s) — L 7 ^^. 
Then the maximizing trajectory passing through an ar- 
bitrary (i, t) can be reconstructed by setting 7(f) = i and 
solving recursively 

7(s-l) = argm fe ax[/i fc ( S -l)-i^ s ( ) 5) ]+^ 7 w( s ) ( n ) 

for s = t,t — 1,...,1. Here argmax^ is the standard 
notation for the value of k that provides maximum to 
the expression in the r.h.s. of (fTTj). 

The algorithmic implementation of the above goes as 
follows. Solve first Eq. © "upstairs" starting from 
given initial conditions and obtain the set of values 
hi(s),fi2(s), . . . ,Hl(s) for all < s < t. Then choose 
a specific point, say (i,t), and restore the path to this 
point going "downstairs," i.e., back in time, by solving 
Eq. (fTl]) step by step. This procedure defines a trajectory 
maximizing the action Aq in Eq. (flU)) . This class of algo- 
rithms is known in the optimization theory as dynamic 
programming, and Eq. © is called the Bellman equation 
(see, e.g., the classical book [13] )■ 

B. BD heaps and the Burgers turbulence 

It turns out that there is a far-reaching analogy be- 
tween the BD deposition model and phenomenology of 
"shocks" and "minimizers" for the Burgers or Hamilton- 
Jacobi equation with random forcing (see, e.g., [26]). We 
first recall the latter. 

Consider the inviscid Burgers equation 

d t u + ud x u = -d x 'n(x 1 1), 

where n(x, t) is the forcing potential. The substitution 
u = d x h transforms this equation into 

d t h+{d x h) 2 /2 + i 1 {x 1 t)=0. 

More generally, one can consider the Hamilton- Jacobi 
equation 

d t h + H(d x h) + n(x,t) = 0, (12) 

where H (p) is a convex function representing the kinetic 
energy. Using the Legendre transform representation 
H(j>) = max„[pw — L(v)], one can write 

d t h + vd x h- L(v) + 77(2;, t) < 




FIG. 2: (a) Heap growing by sequential deposition with highlighted crevices; (b) the growing heap in the (2 + 1) dimensional 
spacetime; (c) density plot of 2nd local difference (discrete analog of 2nd derivative) of the height, which highlights the 
discontinuities corresponding to shocks. Panes (a) and (c) represent front and top views, respectively, of the three-dimensional 
structure in pane (b). 



with equality only for d x h = L'(v), i.e., v — H'(d x h). 
Hence along any trajectory 7(f) the rate of change of h 
is bounded by the Lagrangian 

-/i(7,t)<L(7)-77( 7 ,t) 

(here 7 = dj/dt), which implies for any 7 passing 
through x at time t that 

h(x,t)<Al[j] = h^(0),0)+[ [L( 7 )-q( 7l s)]dfl (13) 

Jo 

with equality only for minimizers of the action, which 
must satisfy the equation 

7 (i) = ff'(0 B / l (7 J i)). (14) 

The Hamilton-Jacobi equation (112[) is thus intimately 
connected with the variational problem of minimizing the 
action (fT3"]) . just as the Bellman equation (O arises in 
maximization of the discrete "action" (fit)]) . Note in par- 
ticular the similar structure of the action (the difference 
in sign results in maximization replacing minimization 
in the discrete case). Moreover, a known solution h to 
(I12p allows to reconstruct minimizing trajectories using 
(Q3J), much as (fTTj) generates maximizing trajectories in 
the discrete problem. 

It is therefore natural to consider the discrete maxi- 
mizing trajectories defined in the previous subsection as 
analogs of continuous minimizers. There is one apparent 
difference: continuous minimizers never cross, while dis- 
crete maximizing paths merge and form tree-like struc- 
tures. However continuous minimizers have a tendency 
to approach each other with exponential rate in reverse 
time due to hyperbolicity, and in the discrete case the 
same hyperbolicity manifests itself in the exponentially 
decreasing probability for two adjacent maximizers to 
stay separate as time runs backwards. 



We are now in position to establish the relation be- 
tween discrete maximizers and connected paths defined 
within the heap in Section [TTJ Lift the maximizing tra- 
jectories to the (i,t,h) space by setting h = Aq(^) for a 
maximizer 7 such that j(t) — i. Then connected paths 
are given by the projection of these "lifted" maximizers 
to the (x,h) plane (see Fig. [2|. In other words, the in- 
tervals of time between successive dropping events along 
a maximizer are collapsed into unit steps in h. Corre- 
spondingly the transversal fluctuations of maximizers as 
a function of time are transformed to transversal fluctu- 
ations of connected path as a function of height h. 

The analogy between continuous minimizers and dis- 
crete maximizers extends to shocks. In the Burgers tur- 
bulence it typically happens that two or more minimizing 
trajectories, which start at different initial locations, pass 
through same point x at time t, so that the map from 
(x y t) to the initial location is discontinuous (see, e.g., 
126]). These discontinuities are called shocks; in space- 
time they form continuous shock curves. This definition 
is obviously parallel to the definition of shocks given in 
the BD setting in Section |H1 (and has inspired the latter). 



IV. FROM BD PATTERNS TO AIRY 
PROCESSES 

A. Basics of classical KPZ scaling 

Recall first the basics of classical KPZ scaling related 
to Airy 1 process, which is closest to our setting. The 
scheme described below is due to Sasamoto [2^ |. 

Consider a directed random walk on a (1 + 1) dimen- 
sional lattice. Suppose that the space-time lattice is 
equipped by a random potential with independent val- 
ues rji(s) at each point (i, s). Then for every i one can 
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consider the maximum of an action over all random walk 
paths of length t terminating at that point, i.e., define 



Oi(t) 



7: l{i)=i 



Ks<t 



V~f{s){s)- 



where Lk.i is a "kinetic" part of the action that ensures a 
certain control of how far the trajectory 7 can jump over 
unit time steps. It is easy to see that i _1 a,(t) — > c at 
t — > 00, where c is some nonrandom constant independent 
of i. We now consider the rescaled process 



A t (x) 



1 



(a at 2/3 x (t) - ct). 



(15) 



The main statement is that A t (x) converges as t — >• 00 
to a universal spatially homogeneous limit process called 
Airyx(x). Universality here means that whenever one 
optimizes in a disordered medium the action of a path 
from a point that varies over a line to a parallel line 
separated from the first one by distance t ( "point-to-line 
last- passage percolation"), the process corresponding to 
the optimal action converges as t — ¥ 00 to the Airyi 
process. 

Note that spatial homogeneity of Airy±(x) immedi- 
ately follows from the construction. Of course one has 
to ensure convergence by subtracting the mean value of 
order t, normalizing the difference by t 1 / 3 and rescaling 
the starting point by t 2 ! 3 . The constants a and /3 in (fTS"]) 
are nonuniversal and should be chosen properly to en- 
sure convergence to the standard Airy process. A simi- 
larly rescaled "point-to-point" percolation results in the 
Airy2 process. 



B. The Airy process for BD pattern 

As we have shown in Section HTT1 the height function in 
the BD process can be viewed as given by maximization 
procedure for random paths in random potential. The 
only difference with the classical picture just described 
is related to the rarity of the deposition events. In other 
words, in order to achieve the displacement of order 1 in 
space direction one needs time of order L. This explains 
why time has to be rescaled. 

The most natural way to do this is through a local 
stochastic change of time variable. Namely we collapse 
the time between two deposition events to 1. This is 
exactly the transformation from (lifted) maximizers to 
connected paths presented in Section [TTTl It is therefore 
no surprise that the Airy\ process can be obtained from 
the BD height function: 

Km ft( f /£)i/3 (K(t/L)V*x(t) - t/L) = Airyxix). (16) 

This formula simply indicates that the appropriately 
rescaled height function in BD is the visualization of the 
process which converges in the thermodynamic limit to 
the Airy process. 



C. The "hairy Airy" process 

The Airyi process carries only part of the informa- 
tion about the system: it is oblivious to the maximizing 
trajectory associated to the (rescaled) point (x,t). It is 
therefore natural to consider the limit 



a atV3 x (t)-ct Jat^xAts) 



fit 1 / 3 



at 2 / 3 



-> (Air yi (x),T x (s)), (17) 



where 7^4 is the maximizing trajectory that passes 
through i at time t and r a; (s) is a continuous path de- 
fined over [0, 1] such that r a (l) = x. We call this limit 
the equipped Airy process. 

As just before, in the BD setting we collapse time in- 
tervals between adjacent deposition events to unit steps 
and get particle paths instead of maximizing trajecto- 
ries in formula (fTTj) above. Applying transversal rescal- 
ing a(t/ L) 2 / 3 x and height rescaling ^(t/L) 1 / 3 as in (|16p . 
we get a realization of equipped Airy process from the 
rescaled BD heap. This process describes the joint distri- 
bution of fluctuations of the height function and transver- 
sal displacements of cluster boundaries in the spatially 
homogeneous BD process. 

In other words, the rescaled height function for the BD 
model alone is a realization of the Airyi process, while 
the rescaled height function together with the rescaled 
forest of maximizers corresponds to a realization of the 
equipped Airy process. The distinctive geometric fea- 
tures of this joint process suggests the name "hairy Airy 
process," cf. Fig. [3^. 

We demonstrate the existence of correlations in the 
joint distribution for the hairy Airy process by comput- 
ing numerically the joint distribution of the height fluc- 
tuation Ah at the top of a shock and the corresponding- 
displacement Ax of the cluster boundary. To be pre- 
cise, we compute the correlation coefficient between the 
fluctuations of the displacement Ax of the right bound- 
ary of a cluster (or equivalently a backbone) and the 
height fluctuation Ah at the top right point of the same 
cluster, see Fig. For convenience we explicitly recall 
here the standard definition of the correlation coefficient 
corr{a, b} between two random variables a and b: 



corrja, b} 



((a-(a))(b-(b)) 
^((a-(a)) 2 )((b-(b)y< 



(18) 



Let the top right particle of some cluster be located at 
time t in column j. Let hj(t) be its height and (h(t)} 
the mean height of the whole surface at time t. De- 
note furthermore hj(t) — (h(t)} by Ahj(t) and the dis- 
placement of the right boundary of the same cluster at 
time t, measured from the position of this boundary at 
t = 0, by Axj(t). We fix a time t, collect for each clus- 
ter the joint information (Ahj(t), Axj(t)), and perform 
averaging over all clusters. Behavior of the correspond- 
ing time-dependent correlation coefficient corr{A/i, Ax}, 
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FIG. 3: (Color online) a) Correlation between the height of 
cluster's boundary and the displacement of the corresponding 
connected path; b) the corresponding correlation coefficient; 
c) the averaged difference between the height of cluster's right 
boundary and the mean height of the BD growing interface. 



where the angle brackets correspond to averaging over 
the sample, is shown in Fig. [5Jd for different time values. 

Strong correlations between the vertical and horizontal 
displacements of cluster boundaries are clearly seen in the 
data. The negative sign of these correlations is due to the 
fact that the height of the top right particle in a typical 
cluster is smaller than the averaged height of the growing 
BD interface. This observation is supported by Fig. [3J;, 
where the averaged difference between the height of the 
cluster right boundary and the averaged height of the 
interface is plotted against time. Clearly this difference 



is always negative and tends to from below as t —> oo. 
One may speculate that growth of the left- and ritghmost 
connected paths in a cluster is slower due to screening 
between neighboring clusters. 

In order to better understand the influence of clus- 
ters on the morphological structure of the growing BD 
surface, we also compute the joint distribution of height 
fluctuations in two columns separated by distance 5 = 3 
in lattice units, as shown in Fig. 2k- Two different situa- 
tions are distinguished: i) two test column belong to the 
same cluster (configuration A), and ii) two test columns 
belong to different clusters, i.e., are separated by a shock 
(configuration B). 
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FIG. 4: (Color online) a) Correlation between the heights 
inside the cluster ("in") and separated by a shock ("out"); b) 
the corresponding correlation coefficients. 



Computing the correlation coefficient corr{Aft,/j, Ah m } 
according to (|18D . we see that correlations between Ah% 
and A/12 inside a cluster are stronger than those across 
a shock between different clusters. 
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V. SCALING ANALYSIS OF BD PATTERNS 



Hence, 



"Thinning" of clusters and wandering of their 
boundaries 



Relying on the connection between shocks and bound- 
aries of clusters, we can directly transfer the scaling ar- 
guments of statistics of shocks developed in [26| to the 
scaling analysis of a growing BD heap and determine the 
values of the scaling exponents a and (3 in the dependen- 
cies c(h) ~ h~ a and Ax 2 (h) ~ hP defined correspond- 
ingly in Eqs (j6|) and ([7]). Recall that c(h) is the averaged 
number of clusters percolating to height h and Ax 2 (h) is 
the mean square displacement of a cluster boundary at 
height h. 

Denote by d(t) the horizontal size of a cluster at time t. 
At t = the cluster has zero size, i.e., d(0) = 0. In 
what follows we shall use the obvious fact that the growth 
time t in the sequential deposition process is proportional 
to the average height h of the growing heap and, conse- 
quently, to the cluster height — see, for example, Fig. [2b- 

The typical value of d(h) can be obtained by scaling 
considerations. Namely, growth of d(h) is determined by 
two additive effects. On the one hand, there is a "driv- 
ing force" promoting the "smearing" of the cluster due 
to the velocity fluctuations. For BD this effect can be es- 
timated as follows. Consider clusters with size of order d. 
Under the uniform random "rainfall" of deposited parti- 
cles, one cluster can randomly screen part of its neigh- 
bors, and increase its own "spot." Since different clusters 
are correlated weakly, it is natural to conjecture that the 
typical scale of fluctuations of cluster sizes is of order of 
Vd. Thus the rate v of cluster "smearing" due to these 
fluctuations is v ~ d/Vd ~ d- 1 / 2 . Speaking more care- 
fully, the above means that the average growth rate of 
the cluster of size d is 



d(t) ~ < 2 / 3 



(20) 



k<j<k+d 



hj), 



where k and k + d are the left and the right boundaries 
of some cluster. The increments of hj are uncorrelated 
for the uniform ballistic "rain" and (hj+i — hj) = 0. It 
is therefore natural to expect that v ~ d~ 1 / 2 as conjec- 
tured. 

On the other hand, there is "smearing" of clusters due 
to the random deposition of new particles near the cluster 
boundary. This process can be interpreted as "diffusion" 
of the boundary. Over time t this diffusion leads to the 
smearing of the cluster's horizontal size on typical scale 
of order of \/t. 

The typical size of a growing cluster at time t is deter- 
mined by additive contributions of these two effects: 



d(t) ~tv + t 1 ' 2 




(19) 



The dominant contribution to d(t) comes from the first 
term, which is consistent with the physical intuition. 



Since t ~ h, we immediately come to the conclusion that 
d(h) ~ /i 2 / 3 . (This estimate is a direct paraphrase of the 
arguments provided in 26] for scaling analysis of statis- 
tics of shocks in the (1 + 1) dimensional Burgers equation 
with random forcing). 

The density c(h) of independent clusters surviving up 
to the height h is inversely proportional to the cluster 
size, c(h) ~ [d(h)}- 1 . Thus, 



c(h) ~ h~ 2 '\ 



(21) 



which gives a = 2/3. 

Furthermore, the typical horizontal mean square dis- 
placement (Ax 2 (h)) of a cluster boundary at height h 
can be estimated simply as 



which gives (3 



(Ax 2 {h)) = d 2 {h) 
- 4/3. 



h A /\ 



(22) 



B. Mass distribution of clusters 

This Section contains the scaling analysis of the prob- 
ability P(m) ~ mT T to find a cluster of mass m in a 
large aggregate. To begin with, note that the number of 
particles, i.e., the "mass" m(h) of a cluster percolating to 
height h can be obtained integrating the horizontal size, 
d(h), of cluster at a given height: 



m(h) 



(h') a dti 



(23) 



From (|2ip we know that the cumulative probability of 
clusters surviving until height h is of the order of h~ a . 
This implies that the probability density of clusters at 
height h scales as dh. In order to calculate 

mass distribution, we change variables from h to m, take 
into account that dm ~ h a dh, or dh ~ h~ a dm = 

m — e»/(a+l) anc J gg-j- 



m _ Q /( Q+ l) dm 



P(m) - m- T dm~/i> +1 >d/i~ 

= m -( 2a + 1 )/( a + 1 ) dm. 

' (24) 
For a — 2/3 we get t = 7/5. The exponent t = 7/5 is 
well supported by our own numerical simulations shown 
in Fig. O and, as mentioned in Introduction, it has been 
found in independent laboratory experiments on clus- 
ter formation in quasi-two-dimensional electrochemically 
formed silver branching structures (23|. 

One can say that he clusters are ranked (ordered) 
according to their masses and m is the corresponding 
rank. Thus Eq. (j2"4"|) has similarity with the Zipf's law 
that appears in many areas of science ranging from word 
statistics in linguistics 29] to nuclear multifragmenta- 
tion [30l . l3lj ] where clusters have power-law distribution 
in sizes (masses, charges etc.). 
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m 

FIG. 5: The probability P(m) to find in a large system a 
cluster of mass m. 

VI. CONCLUSION 

In this paper we analyze the internal structure of the 
heap formed in the course of standard homogeneous bal- 
listic deposition with next-nearest-neighboring (NNN) 
interactions in a box. We have paid the most attention to 
the statistics of clusters and the channels (crevices) sepa- 
rating them. We have demonstrated that the BD process 
can be naturally described in terms of "dynamic pro- 
gramming" language associated with the so-called Bell- 
man equation. The "dynamic programming" point of 
view allows systematic translation of the study of clus- 
ters and crevices in the NNN ballistic deposition into 
the language of maximizers and shocks in discrete equa- 
tions of the Burgers or Hamilton-Jacobi type. This is 
the key point of our work. A detailed examination of the 
corresponding continuous limit will be the subject of a 
forthcoming publication. 

In particular, the results of the work (2(| concerning 
the statistics of shocks in (1 + 1) dimensional Burgers 
turbulence with random forcing allow the direct inter- 
pretation for statistics of cluster's boundaries (crevices) 
of growing heap. This connection between shocks and 
crevices has permitted us to compute the scaling expo- 
nents a (a = 2/3) in the dependence (c(h)) ~ h~ a , 
where (c(h)} is the average number of clusters, surviv- 
ing up to height h and /3 (/3 = 4/3) for the mean square 
displacement (Ax 2 (h)) ~ h@ of crevices as a function of 
the height h of the heap. 

We have also extended the scaling analysis to the com- 
putation of the critical exponent r (r = 7/5) of the mass 
distribution of clusters, P(m) ~ m~ T , where P(m) is 
the probability density of clusters of mass m (see Fig. [5] 



for comparison of numerical simulation with scaling de- 
pendence (HU)). The exponent r coincides with the one 
found in real experiments on cluster formation in quasi- 
two dimensional electrochemically formed silver branch- 
ing structures (23| . 

The investigation of the morphological structure of sur- 
face of the growing heap splitted in clusters, has lead 
us to the definition of a new "equipped" Airy process 
for BD, named the "hairy Airy process." In our pre- 
liminary investigation we have analyzed numerically its 
two-point correlation function and have shown the exis- 
tence of essential correlations between the fluctuations of 
the displacement Ax of the cluster's left boundary and 
the height's fluctuation Ah in the top point of the same 
cluster's left boundary location, see Fig. [3^. 

We believe that the described connection between 
crevices and shocks could be a useful tool for deeper 
understanding of both topics, NNN ballistic deposition 
and Burgers turbulence. For NNN ballistic growth we 
could apply the machinery developed in turbulence, while 
for turbulence we could use NNN ballistic deposition for 
straightforward visualization of some complex chaotic be- 
havior. 

Let us end up by noting that many important and puz- 
zling questions concerning the growth of the heap have 
not been touched in this paper. For instance, we have 
not discussed the question mentioned in the Introduc- 
tion: why the bulk density of the heap does not coincide 
with the density of local maxima of the growing surface. 
Our guess is that the discrepancy between these densi- 
ties, Pbuik and Psurf is due to the presence of crevices in 
the heap. Another example remaining almost without 
the attention deals with the consideration of aging in the 
growing heap. The investigation of the correlation be- 
tween two heights inside a cluster and separated by the 
crevice considered in Fig. 0] gives some hint about the 
aging of the heap, however we have not considered the 
correlation between two heights separated by the crevice 
of finite depth. 
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